﻿#include "MeshTool.h"
#include "Algorithm.h"
#include "Mesh.h"
#include "TopoTables.h"
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
//#include <vtkm/worklet/Triangulate.h>
#include <vtkm/filter/geometry_refinement/worklet/Triangulate.h>
#include <vtkm/worklet/WorkletCellNeighborhood.h>

namespace MeshTool
{
namespace detail
{
struct NumFacesPerCell : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn cell_set, FieldOut numFacesInCell);
  using ExecutionSignature = void(CellShape, _2);

  template<typename CellShapeTag>
  VTKM_EXEC void operator()(CellShapeTag shape, IdComponent& num_faces) const
  {
    num_faces = GetNumberOfSides(shape);
  }
};

struct BuildDuplicateFace : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn cellset, FieldOut face_id, FieldOut face_info);
  using ExecutionSignature = void(CellShape, PointIndices, InputIndex, VisitIndex, _2, _3);
  using ScatterType = ScatterCounting;

  template<typename CellShapeTag, typename C2NType, typename FaceInfo>
  VTKM_EXEC void operator()(CellShapeTag cell_shape,
                            const C2NType& c2n,
                            Id inputIndex,
                            IdComponent visitIndex,
                            Id3& faceId,
                            FaceInfo& face) const
  {
    CellFaceCanonicalId(visitIndex, cell_shape, c2n, faceId);
    face = { inputIndex, static_cast<UInt8>(visitIndex) };
  }
};

struct BuildUniqueFace : public WorkletReduceByKey
{
  using ControlSignature = void(KeysIn faceKey,
                                ValuesIn faceInfo,
                                //ReducedValuesOut face_link,
                                WholeArrayOut cell_links);
  using ExecutionSignature = void(_2, /* _3,*/ _3);

  template<typename FaceInfoPair, typename CellToCellType /*, typename FaceInfo*/>
  VTKM_EXEC void operator()(const FaceInfoPair& face_info,
                            //Vec<FaceInfo, 2>& face_link,
                            CellToCellType& cell_link) const
  {
    if (face_info.GetNumberOfComponents() == 2) // 内部面
    {
      //face_link[0] = face_info[0];
      //face_link[1] = face_info[1];

      auto cell_id = face_info[0].Get().first;
      auto face_id = face_info[0].Get().second;
      cell_link.Get(cell_id)[face_id] = face_info[1];

      cell_id = face_info[1].Get().first;
      face_id = face_info[1].Get().second;
      cell_link.Get(cell_id)[face_id] = face_info[0];
    }
    else if (face_info.GetNumberOfComponents() == 1) // 边界面
    {
      //face_link[0] = face_info[0];
      // TODO: 引入边界条件
      //face_link[1] = { -1, 0 };

      auto cell_id = face_info[0].Get().first;
      auto face_id = face_info[0].Get().second;
      cell_link.Get(cell_id)[face_id] = { -1, 0 };
    }
    else
    {
      this->RaiseError("face_duplicate.GetNumberOfComponents error.");
    }
  }
};

struct ComputeTopolgy : WorkletVisitCellsWithPoints
{
  using ControlSignature = void(CellSetIn, FieldInPoint, FieldOut, FieldOut);
  using ExecutionSignature = void(_2, _3, _4);

  template<typename PType>
  VTKM_EXEC void operator()(const PType& pts, vtkm::Vec3f& centroid, vtkm::Vec3f& spacing) const
  {
    if (pts.GetNumberOfComponents() == 4)
    {
      centroid = (pts[0] + pts[1] + pts[2] + pts[3]) / 4;
      spacing[0] = (pts[1] - pts[0])[0];
      spacing[1] = (pts[2] - pts[1])[1];
      spacing[2] = 1;
    }
    else if (pts.GetNumberOfComponents() == 8)
    {
      centroid = (pts[0] + pts[1] + pts[2] + pts[3] + pts[4] + pts[5] + pts[6] + pts[7]) / 8;
      spacing[0] = (pts[1] - pts[0])[0];
      spacing[1] = (pts[2] - pts[1])[1];
      spacing[2] = (pts[4] - pts[0])[2];
    }
    else
    {
      RaiseError("unfinished");
    }
  }
};

struct BuildBCMape : vtkm::worklet::WorkletCellNeighborhood
{
  using ControlSignature = void(CellSetIn, FieldOut bnd_id);
  using ExecutionSignature = void(Boundary, _2);

  template<typename BndIdType>
  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary, BndIdType& bnd_id) const
  {
    bnd_id = { -1, -1, -1, -1 };
    if (boundary.IJK[0] == 0) // left, inlet
    {
      bnd_id[3] = 3;
    }
    if (boundary.IJK[0] == boundary.PointDimensions[0] - 1) // right, outlet
    {
      bnd_id[1] = 1;
    }
    if (boundary.IJK[1] == boundary.PointDimensions[1] - 1) // top, wall
    {
      bnd_id[2] = 2;
    }
    if (boundary.IJK[1] == 0) // bottom
    {
      bnd_id[0] = 0;
    }
  }
};

struct BuildBCType : vtkm::worklet::WorkletMapField
{
  using ControlSignature = void(FieldIn element_side, FieldIn bc, WholeArrayOut cell_links);
  using ExecutionSignature = void(_1, _2, _3);

  template<typename FaceInfoPair, typename CellToCellType>
  VTKM_EXEC void operator()(const FaceInfoPair& face_info,
                            BCType bc_type,
                            CellToCellType& cell_link) const
  {
    auto cell_id = face_info.first;
    auto face_id = face_info.second;
    auto element_side = cell_link.Get(cell_id);
    //std::cout << int(face_id) << std::endl;
    //std::cout << "N: " << element_side.GetNumberOfComponents() << std::endl;
    element_side[face_id] = { -1, static_cast<SideId>(bc_type) };
    //std::cout << element_side[face_id].Get().first << std::endl;
    //[face_id];
    //std::cout << element_side[face_id].Get().first << std::endl;
    //std::cout << /*static_cast<int>*/(element_side[face_id].Get().second) << std::endl;

    //std::cout << cell_id << std::endl;
    //std::cout << face_id << std::endl;
  }
};

}




void AttachBoundaryCondition(Mesh& mesh, const std::map<std::string, BCType>& bnd_name_2_bc_type)
{
  using FaceInfo = typename Mesh::FaceInfo;

  const auto& boundary_info = mesh.GetBoundaryInfo();
  auto& ss_id_to_name = boundary_info.GetIdToName();

  std::vector<FaceInfo> element_side;
  std::vector<BCType> bc_type;

  auto sideset_map = [&ss_id_to_name]()
  {
    std::stringstream ss;
    ss << "\n";
    for (auto& var : ss_id_to_name)
    {
      ss << "                   (" << var.first << " --> " << var.second << ")\n";
    }
    return ss.str();
  };


  console::Info("  Side Set: ID --> Name", sideset_map());

  const auto& side_set = boundary_info.GetSideSet();
  for (auto& val : side_set)
  {
    // 1.
    element_side.emplace_back(val.first, val.second.first);

    // 2.
    auto bnd_id = val.second.second;
    if (!ss_id_to_name.count(bnd_id))
    {
      console::Error("不存在Boundary ID: ", int(bnd_id), "\n", sideset_map());
    }

    auto bnd_name = ss_id_to_name.at(bnd_id);

    if (!bnd_name_2_bc_type.count(bnd_name))
    {
      console::Info("可用的Bounary Name", sideset_map());
      console::Error("不存在Boundary Name: ", bnd_name);
    }
    bc_type.emplace_back(bnd_name_2_bc_type.at(bnd_name));
  }

  ArrayHandle<FaceInfo> element_side_array =
    vtkm::cont::make_ArrayHandle(element_side, vtkm::CopyFlag::On);
  ArrayHandle<BCType> bc_array = vtkm::cont::make_ArrayHandle(bc_type, vtkm::CopyFlag::On);

  std::cout << bc_array.GetNumberOfValues() << std::endl;
  std::cout << element_side_array.GetNumberOfValues() << std::endl;
  Invoker{}(detail::BuildBCType{}, element_side_array, bc_array, mesh.GetCellLinks());
}

void BuildCellLinks(Mesh& mesh)
{
  using FaceInfo = typename Mesh::FaceInfo;
  ArrayHandle<IdComponent> num_face_per_cell;
  Invoker{}(detail::NumFacesPerCell{}, mesh.GetCellSet(), num_face_per_cell);

  ArrayHandle<Id3> face_id;
  ArrayHandle<FaceInfo> face_info; //包含重复face

  Invoker{}(detail::BuildDuplicateFace{},
            ScatterCounting{ num_face_per_cell },
            mesh.GetCellSet(),
            face_id,
            face_info);

  auto& face_num_offset = mesh.GetFaceNumOffset();
  mesh.GetFlatCellLinks().Allocate(face_info.GetNumberOfValues());
  face_num_offset = vtkm::cont::ConvertNumComponentsToOffsets(num_face_per_cell);
  num_face_per_cell.ReleaseResourcesExecution();

  auto cell_links = mesh.GetCellLinks();
  Keys<Id3> face_keys(face_id);
  face_id.ReleaseResourcesExecution();
  Invoker{}(detail::BuildUniqueFace{}, face_keys, face_info, cell_links);
  face_info.ReleaseResourcesExecution();
}

vtkm::cont::DataSet CreateStructuredMesh(const Id2& point_dim,
                                         const Vec2f& orign,
                                         const Vec2f& space)
{
  vtkm::cont::DataSet data_set = vtkm::cont::DataSetBuilderUniform::Create(point_dim, orign, space);

  ArrayHandle<Vec3f> spacing, center;
  //ArrayHandle<Id> b_id;
  //ArrayHandle<vtkm::Int16> blk_id;

  vtkm::cont::Invoker{}(detail::ComputeTopolgy{},
                        data_set.GetCellSet(),
                        data_set.GetCoordinateSystem(),
                        center,
                        spacing);

  data_set.AddCellField("spacing", spacing);
  data_set.AddCellField("center", center);
  //data_set.AddPointField("boundary_id", b_id);
  //data_set.AddCellField("blk_id", blk_id);

  return data_set;
}

vtkm::cont::DataSet CreateStructuredMesh(const Id3& point_dim,
                                         const Vec3f& orign,
                                         const Vec3f& space)
{
  vtkm::cont::DataSet data_set = vtkm::cont::DataSetBuilderUniform::Create(point_dim, orign, space);

  ArrayHandle<Vec3f> spacing, center;
  ArrayHandle<Id> b_id;
  //ArrayHandle<vtkm::Int16> blk_id;

  vtkm::cont::Invoker{}(detail::ComputeTopolgy{},
                        data_set.GetCellSet(),
                        data_set.GetCoordinateSystem(),
                        center,
                        spacing);


  data_set.AddCellField("spacing", spacing);
  data_set.AddCellField("center", center);

  return data_set;
}


void TriangulateMesh(Mesh& mesh)
{
  vtkm::worklet::Triangulate triangulate;

  auto out_cellset = triangulate.Run(mesh.GetCellSet());
  mesh.GetCellSet() = out_cellset;
  //this->_cell_set = output.GetCellSet();
  //ArrayHandle<Vec3f> pts;
  //Invoker{}(CoodConvert{}, output.GetCoordinateSystem(), pts);
  //this->_coord = CoordinateSystem("coord", pts);
}

}


#include <libmesh/boundary_info.h>
#include <libmesh/elem.h>
#include <libmesh/mesh.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>

namespace MeshTool
{

void Build(Mesh& my_mesh, libMesh::MeshBase& mesh)
{

  mesh.print_info();



  using CoordType = vtkm::Vec3f;
  std::vector<CoordType> coords;
  for (const auto& node : mesh.node_ptr_range())
  {
    //std::cout << *node << std::endl;
    coords.emplace_back((*node)(0), (*node)(1), (*node)(2));
  }

  std::map<libMesh::ElemType, vtkm::CellShapeIdEnum> libmesh2vtkm{
    { libMesh::ElemType::EDGE2, vtkm::CellShapeIdEnum::CELL_SHAPE_LINE },
    { libMesh::ElemType::TRI3, vtkm::CellShapeIdEnum::CELL_SHAPE_TRIANGLE },
    { libMesh::ElemType::QUAD4, vtkm::CellShapeIdEnum::CELL_SHAPE_QUAD },
    { libMesh::ElemType::TET4, vtkm::CellShapeIdEnum::CELL_SHAPE_TETRA },
    { libMesh::ElemType::HEX8, vtkm::CellShapeIdEnum::CELL_SHAPE_HEXAHEDRON },
    { libMesh::ElemType::PRISM6, vtkm::CellShapeIdEnum::CELL_SHAPE_WEDGE },
    { libMesh::ElemType::PYRAMID5, vtkm::CellShapeIdEnum::CELL_SHAPE_PYRAMID },
  };

  std::vector<vtkm::UInt8> shapes;
  std::vector<IdComponent> numindices;
  std::vector<Id> conn;
  for (const auto& elem : mesh.element_ptr_range())
  {
    try
    {
      shapes.push_back(libmesh2vtkm.at(elem->type()));
    }
    catch (const std::exception& e)
    {
      console::Error("不支持的单元类型", elem->type());
    }

    numindices.push_back(elem->n_nodes());
    for (size_t i = 0; i < elem->n_nodes(); i++)
    {
      conn.push_back(elem->node_id(i));
    }
  }
  vtkm::cont::DataSetBuilderExplicit dsb;
  my_mesh.GetDataSet() = dsb.Create(coords, shapes, numindices, conn, "coordinates");

  auto& boundary_info = mesh.get_boundary_info();


  auto& _boundary_info = my_mesh.GetBoundaryInfo();
  for (auto& side_set : boundary_info.get_sideset_map())
  {
    _boundary_info.AddSide(side_set.first->id(), side_set.second.first, side_set.second.second);
  }

  for (auto& name_map : boundary_info.get_sideset_name_map())
  {
    _boundary_info.SideSetName(name_map.first) = name_map.second;
  }
}

}
